Hopping-resolved electron-phonon coupling in bilayer graphene 
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In this paper we investigate the electron-phonon coupUng in bilayer graphene, as a paradigmatic 
case for multilayer graphenes where interlayer hoppings are relevant. Using a frozen-phonon ap- 
proach within the context of Density Functional Theory (DFT) and using different optical phonon 
displacements we are able to evaluate quantitatively the electron-phonon coupling m associated with 
each hopping term 7i. This analysis also reveals a simple scaling law between the hopping terms 7^ 
and the electron-phonon coupling m which goes beyond the specific DFT technique employed. 



I. INTRODUCTION 

Since its discovery, a formidable amount of work has 
been devoted to investigate the electronic and struc- 
tural properties of single-layer and multi-layer graphenes. 
The electron-phonon interaction has in particular at- 
tracted a huge interest for its role in controlling the 
charge transport^"— for providing a powerful interface 
between electronic and structural properties, and also 
because phonon resonances in Raman and infrared spec- 
troscopies, triggered by the electron-phonon interaction, 
represent a useful way to characterize graphenic samples 
and to reveal interesting unconventional effects j^"— 

On the theoretical level, tight-binding (TB) models are 
of fundamental importance since they have been shown 
to catch almost all the electronic features in this systems. 
The vast majority of works based on TB in literature^ 
employ a simple two-parameter TB model, where only 
the nearest neighbor in-plane hopping 70 and the nearest 
neighbor vertical hopping 71 are considered, although, 
when needed, higher order TB terms are included to re- 
produce more detailed features. ^^"^^ Most important in 
bilayer graphene, borrowing the terminology from bulk 
graphite)^ on the atoms Bl and A2, and the hopping 
terms 73, operative between the atoms A1-B2, and 74 
between atom couples A1-A2 and B1-B2 (see Fig. [1]). 
Such terms represent thus the basilar ingredients to build 
a TB model in multilayer graphene with both Bernal 
(ABAB. . . ) and rhombohedral (ABCA. . . ) stacking. 
For instance, the hopping term 73 was shown to be re- 
lated to the trigonal warping and,^^ in bilayer systems, to 
the generation of new Dirac points at finite momentum 
close to the K point. 

Tight-binding models are also widely employed to in- 
vestigate the electron-phonon interaction. Focusing on 
the single-layer graphene, the most relevant in-plane lat- 
tice vibrations are related to the modulation of the near- 
est neighbor hopping 70.— Within this context, for in- 
stance, the optical properties on the E2g phonon band 
at w « 0.2 eV have been throughout investigated^^"'^^ 
as well as the effects on the electronic structure of 
the long wavelength acoustic modes associated with the 



ripples.-^iS, Modeling of the electron-phonon interaction 
in multilayer graphene is also commonly discussed on the 
basis of the modulation of the 70 hopping. Among other 
things, this kind of analysis was useful to show the ro- 
bustness of the Dirac points^^ upon lattice distortions, in 
single- as well in multi-layer systemsi^i^ 

Alternative to tight-binding model. Density Functional 
Theory (DFT) calculations permit to include all the dif- 
ferent kinetic (e.g. hopping) terms at the same level. 
It also permits to provide a quantitative estimate of the 
electron-phonon coupling. Pivot in this context is the 
concept of deformation potential, i.e. the shift of the 
electronic levels upon a frozen phonon lattice distortion, 
which is strictly related to the magnitude of the electron- 
phonon interaction.'^^ In the context of graphenic mate- 
rials, frozen phonon DFT calculations were employed to 
quantify in single layer systems the electron-phonon cou- 
pling associated to the modulation of the 70 term upon 
a lattice distortion u^^^ It can be shown indeed that 
the in-plane optical mode £'29 induces a linear splitting 
Ae of the Dirac states at K point, Ae oc 6au, where a 
is related to the linear coupling of the electronic states 
with the q = E2g mode.— In a TB model, defin- 
ing u as displacement per atom, one gets a = d-ja/du. 
DFT calculations in single-layer graphene obtain a = 4.5 
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FIG. 1: (color online) (a) Atomic structure of multilayer 
graphene with Bernal stacking showing the relevant hopping 
terms ji. Atoms Bl and A2, connected by vertical 71, denoted 
by darker colors, contain also a local crystal field potential. 



eV/Aj^ A similar value a « 4.4 eV/A is found also in 
graphite^Si^ where the linear energy splitting Ae at the 
H point (where the interlayer hopping is unaffective) can 
be shown to be uniquely related to a = d-jo/du. In both 
cases, in graphene and graphite, a GW theory leads to 
slight larger values a — 5.1 ~ 5.3 eV/Ai^ 

Despite large effort has been devoted thus in litera- 
ture to study the electron-phonon interaction related to 
the 7o term, virtually no work has been addressed so far 
to provide a quantitative estimate of the electron-phonon 
coupling associated with the modulation of the other hop- 
ping terms. A quantitative insight on this issue, on the 
other hand, becomes increasing important because of the 
role of such terms to many effects, from the establishment 
of unconventional anisotropic phases in strained bilayer 
systems- to the evaluation of the optical properties of 
the in-plane and out-of-plane phonon mode in multilayer 
systems and in graphite."*^ 

Aim of the present paper is to fill this gap and to pro- 
vide, with a first-principle DFT calculation, a quanti- 
tative study of the electron-phonon coupling associated 
with the modulation of other main hopping terms, both 
for in-plane and for the out-of-plane vibrations. We ad- 
dress this issue focusing on the optical phonon modes at 
q = in bilayer graphene. The modulation of each hop- 
ping term with the relative distance, however, permits 
to provide a generalization of the present results at any 
finite q. 



II. FROZEN PHONON ANALYSIS 
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In this paper we consider single-layer and bilayer 
graphene with typical Bernal stacking. We take the in- 
plane nearest-neighbor carbon-carbon distance b = 1.42 
A (a = 2.46 A the lattice constant), and the inter- 
layer distance d = 3.35 A. Such lengths rule thus the 
magnitude of the in-plane hopping term 70 and out-of- 
plane hopping terms 7^ on the relative distance of the 
corresponding atoms. For sake of simplicity, we denote 
with bi the distance associated with each hopping term 
in the undistorted structure, namely 69 — |Rai ~ R-Bil, 
61 = |Rbi -Ra2|, ^3 = |Rai -Rb2|, 64 = |Rai -Ra2|. 
We assume that on the local scale the hopping terms 7^ 
depend uniquely on the relative distance r, 7^ = 7i(?')- 
The modulation of such hopping terms induced by the 
lattice displacement determines thus the electron-phonon 
interaction. In full generality, we define thus a electron- 
phonon coupling as ai = — c?|7i|/c?r|r=6i Note that, since 
the amplitude of the hopping parameters |7i(r)| generally 
decreases with increasing the distance r, we have intro- 
duce an explicit sign (-) in the definition of ai so that the 
corresponding electron-phonon couplings is by definition 
chosen to be positive. 

In order to reveal the electron-phonon coupling ai for 
each hopping parameter 7^, we consider the i?2g mode 
for the single layer graphene, and the -Big^ , ^'292 ^.nd 
E2gi , for the bilayer graphene as sketched in Fig. [5^. 



FIG. 2: (color online) (a) Sketch of atomic displacements for 
the relevant phonon lattice modes here considered in single- 
layer graphene (£25), and in bilayer graphene (-B191, i?2g2i 
£-291, Eiu)- (b) Representative band structure of single layer 
and bilayer graphene upon a frozen phonon lattice distortion 
(red dashed lines) with it = 0.1a = 0.0246 A for the E2g mode 
in single layer graphene (left panel) and for the £^292 mode in 
bilayer graphene (right panel). Also shown is the undistorted 
band structure and the band labels (ei-e4 from top to bottom 
band. 



We focus on the deformation potential close to the K 
point, where one-particle low-energy excitations are in- 
volved, which makes a DFT approach particularly efh- 
cient. We compute the electronic band structure in the 
presence of a static frozen phonon displacement by us- 
ing a plane-wave implementation^^ of the density func- 
tional theory in the local-density approximation (LDA) 
for the exchange-correlation potential. Ultra-soft pseu- 
dopotential for carbon was used with plane- wave (charge 
density) cutoff of 40 (400) Ry. A uniform wave-vector 
grid of 18x 18 in the irreducible Brillouin-zone with cold- 
smearing of 0.02 Ry was sufficient to converge the calcu- 
lated quantities to the required accuracy. 

In order to provide a common framework for all the lat- 
tice modes of single-layer graphene and well as of bilayer 
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graphene/graphite, we analyze the deformation poten- 
tial as function of m, where u represents the magnitude 
of the lattice displacement of each atom. We consider 
both degenerate in-plane modes along x and y directions, 
and the non degenerate out-of-plane modes. For each 
case we choose, respectively, uai,x — Ua;X, uai,j^ — Uyy, 
uai,z = UzZ. The displacement of the other atoms is thus 
univocally determined by the components of the wavevec- 
tor of the phonon mode. 

Representative electronic structures of the single-layer 
and bilayer graphene in the presence of lattice distor- 
tions are shown in Fig. ^jp. Focusing at the K point we 
can expect, according the different modes considered, an 
opening of a gap for the Dirac energy levels and a fur- 
ther modulation of the upper and lower energy bands. 
In the bilayer system, we label the four 7r-bands as £1-64, 
from the top to bottom energy, as shown in Fig. [2}d, and 
we denote Ae the possible splitting of the Dirac state 
Ae = £2 — €3 and = £1 — £4 the energy difference be- 
tween the upper and lower band. We fix for convenience 
the energy zero of our band structure at the Dirac point 
of the undistorted system. It is important to stress that 
our procedure indeed involves only energy differences so 
that the absolute energy position of the band structure 
is irrelevant. 

It is also useful to introduce here the low-energy Hamil- 
tonian for the undistorted lattice structures. Using stan- 
dard notations, the single-layer and bilayer graphenes are 
thus described respectively by the Hamiltonians: 
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where /k = e~'''=^''/^-|-2e*''"°/^^cos(fcj^a/2), and where 
8 is the difference of the crystal field probed by the Bl- 
A2 carbon atoms in the bilayer structure with respect to 
the A1-B2 atoms. 

The band structure for the undistorted bilayer 
graphene is shown in Fig. ^p. Equating the TB ana- 
lytical expressions with the computed DFT eigenvalues 
we get £i = (5 + 71 = 0.3620 eV, £4 = ,5 - 71 = -0.3382 
eV, which permits to evaluate the parameters 6 = 0.0119 
eV and 71 = 0.3501 eV. 



A. Single-layer graphene 
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FIG. 3: Splitting of the Dirac point Ae as a function of the 
E2g lattice distortion (u) in single layer graphene as evaluated 
from frozen phonon DFT calculations. The slope of Ae vs. u 
gives Qfo = 4.41 eV/A. 



x-direction we have thus: 

r-B2„, X / 



if--, 



70 /k 



70 /k + aog^u,^ 







,(3) 



where 5k = 2e-'''-''/^ - 2e'''''''/'^^ cos{kya/2). At the 
K— 47r/3a(0, 1) point, we get /k = 0, = 3, so that 



Ao,: 
Ao,:. 



(4) 



where Aq^^ = SaoUx- The degenerate levels at the 
Dirac point result thus splitted in single layer graphene 
upon a E2g lattice distortion along the x-axis of a to- 
tal amount Ae^^!>{u) = 2|Ao,a;| = 6q!o| Ux\, in agreement 
with Refs^>^^. 

A similar result can be obtained by considering lattice 
displacements along the y-direction. In this case we have: 





7o/k + aohluy 



'k ^ "0"k"y 

where = i2V3e*'=-"/2^sin(A:j,a/2). 
47r/3a(0, 1) point, we get Hk — 3i, so that 



70 /k + aohwUy 




,(5) 



At the K= 





-jAo,y 



iAo.y 




(6) 



with Ae^^!>{u) — 2|Ao,j,| = 6ao|wj,| also in this case, re- 
flecting the perfect degeneracy of the x vs y in-plane lat- 
tice vibrations. 

Our DFT computed splitting is shown in Fig. [3l from 
which we get ao = 4.41 eV/A, in nice agreement with 
Refs. IsgHssI . We found virtually no difference for lattice 
distortions along x ot y direction on this range of u. 



1. E2g mode 



With these notations, we can now consider, as a pre- 
liminary check, the frozen phonon Hamiltonian of the 
single-layer graphene upon the E2g distortion. Along the 



B. Bilayer graphene 

Once evaluated the in-plane electron-phonon coupling 
ao associated with the 70 hopping term in the single- 
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FIG. 4: (a) Energy diflFerence AS^^'i (u) between upper and 
lower bands as a function of the Big^ lattice displacement itz 
in bilayer graphene. (b) Corresponding ratio AB^^fi /u as 
a function of the lattice distortion u. The extrapolation for 
u->-0 gives ai = 0.608 eV/A. 



distortion. We can note that, as a consequence of the 
symmetry preserving displacements, no gap is opened 
at the K point. Useful information is however encoded 
in the energy difference between the high-energy bands 
E^^ni [u) — ei(u) — e4,{u) which, from Hamiltoniar(51 re- 
sults 



E'^isi (u) = 271 + 4aiw. 



(10) 



We can evaluate thus the electron-phonon coupling ai 
from the linear dependence of /S.E^^n (u) = [E^^n (u) — 
E^^si (0)] on u. The calculated DFT dependence of 
AE'^isi (u) as a function of the vertical displacement Uz 
is shown in Fig. S^, whereas the ratio IS.E^^o'i- {u) /u is 
shown in Fig. Hb, whose extrapolation for — gives 
ai ^ 0.608 eV/A. 



layer graphene, we can now address the role of higher 
order hopping terms in multilayer graphenes, using the 
bilayer graphene as a suitable tool. 

1. Big^ mode 

We first consider the out-of-plane Big^ mode, as de- 
picted in Fig. [2^. This is a quite peculiar mode since it 
does not lift any symmetry of the crystal. We can thus 
still write the four energy levels at the K point as 

£2/3 = 0, (7) 
£1 = S{u)+j,{u), (8) 
£4 = S{u)-j,{u), (9) 

where we have explicitly expressed the intrinsic depen- 
dence of the parameters S and 71 on the lattice 



2. E2g2 mode 

The mode is quite peculiar as, since it does not lift 
any symmetry of the system, it does not split the Dirac 
energy levels at the K point. We have shown above how- 
ever that the splitting of the additional upper and lower 
bands can be used to estimate the electron-phonon cou- 
pling associated with 71. Things are richer when other 
modes, reducing the symmetry of the crystal, are con- 
sidered. In this case useful information about different 
electron-phonon coupling are encoded in the splitting of 
the Dirac point as well as in the u-dependence of the 
differences between high-energy bands, AE{u). 

Let us consider for instance the electronic structure of 
the bilayer graphene under a i?2g2 lattice distortion. If 
we consider only the leading order linear electron-phonon 
couplings ai, we can thus write 
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70 /k + a^glu^ 
\ 73/k - cosOasgkUx 
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73/k ~ cos Ba^g^u^ 
74/k 
70 /k + aogkUx 




(11) 



where cos 9 = b/VlP + cP « 0.39 is a geometric factor 
accounting for the projection of the lattice displacement 
along the direction of the 73 bond. Evaluated at the K 
point, it reads 
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where Ag^^/^ = 3cos9a3Uj;/y. 



The eigenvalues e's can be thus obtained from the sec- 
ular equation: 

e^iS - ef + 2A2 - e) - e'jf + 2Ag^,A3,.7i 
+Al,-Al,i6-ef + Al,jf = 0, (13) 

Eq. predicts a linear splitting of the Dirac levels as 
a function of Ux- Linearizing with respect to Ux we find: 

Ae^^»= = 2|A3,,| = 6cos0a3|u,|, (14) 
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FIG. 5: (a) Splitting of the Dirac point Ae as a function of 
the ^^292 lattice distortion u in bilayer graphene. Empty cir- 
cles are DFT frozen phonon calculations for u^, while empty 
squares for Uy. Solid lines are the solution of the correspond- 
ing analytical models in Eqs. p3p . (|16p . using 5 — 0.0119 eV, 
71 = 0.35 eV and a-j = 0.54 eV/A. The value of qq eV/A is 
irrelevant on this quantity in this range, (b) Corresponding 
Ae/|M| as a function of the lattice distortion u. The extrapo- 
lation of the DFT data for m — > gives an unbiased estimate 
of Q3 = 0.54 eV/A. 



which permits us to evaluate from the linear splitting 
at the K point of the Dirac bands in bilayer graphene 
upon a i?2g2 lattice distortion. In Fig. [3^ we show 
the linear splitting Ae^^"^ computed by using our frozen 
phonon DFT calculations for different (open circles). 

The linear extrapolation of Ae-^^"^ /u^ for dux 0, as 
shown in Fig. [5]d, gives us thus an unbiased estimate of 
aa = 0.54 eV/A. 

The accuracy of such estimate, as well as of the tight- 
binding analysis here considered, can be checked by using 
this last value {as = 0.54 eV/A) and the TB parame- 
ters previously evaluated in an independent way in the 
undistorted structure {6 — 0.0119 eV, 71 — 0.35 eV) to 
calculate the splitting on a wider range of Ux , without the 
linearization, but solving Eq. (fT3l) . The analytical results 
obtained in this way are in excellent agreement with DFT 
calculations proving thus the full intrinsic consistency of 
the value of as with respect the other TB parameters. 

Note also that the DFT calculations predict a critical 
value Ux where the gap at the K point close, reconstruct- 
ing there thus, for this particular value of u^;, a Dirac 
cone. This peculiar feature can also be understood using 
the TB model. As a matter of fact, from an inspection 
of Eq. (fT3|) . one can find two very close critical values 
Ux = —as cos6'(7i ±(S)/3aQ where the gap at the K point 
closes. These points are however so close that they can- 
not be resolved on the scale of Fig. \E\ The reconstruction 
of the Dirac cone at the K point is a mixed combination 
of the effects of the trigonal warping induced by 73 and of 
the additional effects related to the lattice distortion. In 
the undistorted structure, indeed, we know that the effect 
of 73 in bilayer systems is to induce satellite Dirac cones 
at finite k in addition to the conventional one at the K 



point. Lattice distortions induce, as well as in single-layer 
graphene, a shift of the main Dirac point away from the 
K point, opening thus there a gap. The satellite Dirac 
points however move as well as functions of the lattice 
distortion. At a certain value, Ux, one of the satellite 
Dirac points is moved again across the K point, and this 
feature is reflected in the closing of the gap in Fig. [S] 
at a finite Ux- The value of Ux agrees also in excellent 
way with the above analytical estimate from the tight- 
binding model. On the other hand, for Uy displacements, 
the Dirac point moves in an orthogonal direction with 
respect to the K point and no reconstruction of Dirac 
cones at K is possible. A more detailed analysis of this 
issue is provided in Appendix IIIIl 

Finally, as a last check of our analysis, we computed 
also the frozen phonon energy splitting for i?2g2 lattice 
displacements along y. DFT calculations are shown in 
Fig. [5] as empty squares. To extract information about 
the electron-phonon coupling, we analyze the Hamilto- 
nian at the K point which reads now: 
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, (15) 



with a secular equation: 

e\S-er+2Al^e{S-e)-e'j', 



A2 (5-e)2 + A2 7?=0. 



(16) 



Note that, unlike the displacements along x [Eq. P^ ]. 
Eq. (|16p is symmetric with respect to Uj^ — > — 



small values of Uy, we once more obtain 
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:-E2g2 _ 



2|A3^y| = 6cos6'a 



3\Uv\ 
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(17) 



reflecting the degeneracy, at the linear level, of the £'232 
mode along the two directions. The extrapolation of 
Ae/uy coincides with Ae/ux for u — )■ 0, providing thus 
the same value as = 0.54 eV/A. 

It is also interesting to give a look now at the depen- 
dence of E^^^ at the K point with respect to the lattice 
displacement Ux- For these levels we find a quadratic de- 
pendence on Ux- Expanding Eq. p6p at the second order 
with respect to Ux, we get 



ei 



S + ji + 



71 + (5' 



so that 



£4 = (5 - 71 



AE^3'2 ^ 271 



A2 

71 - (5' 

271 At 

7? - 



(18) 
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The DFT calculations (open symbols) of the it- 
dependence of AE'^3^ are shown in Fig. [5] (panel a), as 
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FIG. 6: (a) Variation of AE as a function of the i?292 lattice 
distortion Ux in bilayer graphene. (b) Corresponding ratio 
AE/u'^. The extrapolation for it^, — > gives qq = 4.40 eV/A. 
Empty circles are DFT data, while solid lines are obtained 
from the TB model with S = 0.0119 eV, 71 = 0.35 eV and 
ao = 4.40 eV/A and ai = 0.68 eV/A. The value of 03 eV/A 
is irrelevant on this quantity in this range. 



well with the ratio AE"^^^ /u^ (panel b). The extrapola- 
tion of Ai?^f2/u^ for Uj; — >■ provides thus an estimate 
ao = 4.40 eV/A which essentially coincides with the value 
extracted in the single-layer graphene. 

It should be here noted that Eq. (PO)) has been derived 
from Eq. (jlip where only the linear terms in u where re- 
tained. Some care is however needed on this regards since 



we are actually investigating here a quadratic dependence 
on u. A careful analysis shows that further corrections at 
the quadratic order in Eq. ([^0]) appear through the ex- 
plicitly dependence of 71 on u. Taking into account the 
geometry of the lattice displacement, one should write 
thus 



271 



271 Ag,, 



7? 



(21) 



The correction coming from ai are however two orders 
of magnitude smaller that the term oc Ao,^; and they are 
here ineffective. 



3. -B291 mode 

After having determined the electron-phonon coupling 
ao, Qfi, as in bilayer graphene from the frozen phonon de- 
pendence of the energy levels at the K point under Big-^ 
and i?2c/2 distortions, we are now aiming to a correspond- 
ing characterization of the last remaining parameter 
associated with the 74 hopping. The most straightfor- 
ward way to probe it, as we are going to see, is to consider 
the i?2gi phonon mode, as depicted in Fig. [21 

Upon distortion along the ^2gi phonon mode, the 
Hamiltonian reads: 



7o/k 74.h + cos 9a4,gku^ Jsfu - cos Oasg^u^ 

H^^'Uu ) = I '^"■^^ S 71 74/k + cos6'Q;45kU:c 1 ^22) 

^ ^ I 74/k + cos 6la45kWx 71 S 7o/k ^ ' 

cosO-ysf]^- asgi^u^ ^Af^+ cosOaigl-Ux 70/k 



Evaluated at the K point, we thus have: 
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(23) 



where A/^^x = icosOuiiUx, and we can write the secular 
equation: 

e^{5 -ef+2 cos^ OAl ^eiS - e) - 6^72 

+2 cos^ dAl x/^^^xli + cos"* eA\ x - cos^ OAl ^iS - ef 

+ cosHAl^^l = 0. (24) 

Eq. p4)l . predicts also, as (fT3)) . a linear splitting of the 
Dirac level upon lattice distortion associated once more 
with tta, i.e. Ae^^"^ = 6 cos6'a3|ua;|. The value of 03 esti- 
mated upon such lattice distortion coincides with the one 



obtained previously using the -^292 mode, corroborating 
thus the analysis. 

More useful information is however encoded in the 
frozen phonon dependence of AE. Such splitting was 
above employed to estimate directly ao from the frozen 
phonon i?2g2 lattice distortion. In the present E2g^ con- 
text, we can see that we still get, although not a direct, 
an indirect estimate of a4 from the u-dependence of AE. 
We can indeed write 

ei = '5 + 71 + ^1^, (25) 

71+0 

and 

AL 

£4 = (5-71-^^, (26) 
71-0 

so that AE is expected once more to presents a quadratic 
dependence on u. Taking into account, just as in the 
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FIG. 7: (a) AE as a function of the i52gi lattice distortion 
Ux in bilayer graphene. Inset: same on a wider u-region. 
Axis labels in the inset are the same as in the main panel, 
(b) Corresponding ratio AE/u^- From the extrapolation for 
Ux — > 0, and taking into account the contribution of the u 
dependence of 71, we can estimate 04 — 0.32 eV/A. Empty 
circles are DFT data, while solid lines are obtained from the 
TB model with S = 0.0119 eV, 71 = 0.35 eV and ai = 0.61 
eV/A and 0:4 — 0.30 eV/A. The value of ao eV/A is irrelevant 
on this quantity in this range. 



4. Other modes (E\u, B\g^, . . .) 

Other optical modes at q = can be in principle con- 
sidered to investigate the deformation potential due the 
electron-phonon interaction. However, they result to be 
not particularly convenient in order to disentangle the 
role of the different electron-phonon couplings associated 
with the different hopping parameters. Once can see for 
instance that the (also shown in Fig. ^ induces a 
quadratic splitting of the Dirac point as a function of m, 
whose curvature depends on the same level on both 
and a4, so that their values cannot be estimated in an 
unbiased way from an extrapolation for m — 0. Similar 
problems appear when considering the splitting of high 
energy bands for Eiu , or the energy splitting (Dirac point 
as well as high-energy bands) for the other modes. Also 
in these cases, the deformation potential results to be a 
mixing of different electron-phonon coupling, making the 
quantitative evaluation of the ai from these modes not 
reliable. We have however checked, on the other hand, 
that the above values estimated from the Sig^ £'2^2 ^''^'^ 
_E2gi modes reproduce the energy differences of the elec- 
tronic bands at the K point upon other different lattice 
modes. 



i?2£(2 case, the quadratic dependence associated with 71, 
we can write thus 



III. DISCUSSION AND CONCLUSIONS 



AS^^i = 271 



(27) 



DFT calculations for AE'^^^ are shown in Fig. [7^ on 
the same u-scale employed for other lattice modes. Due 
to the smallness of such it-dependence, numerical noise is 
here much larger than in previous analyses. A negative 
quadratic curvature can be however still clearly observed, 
which is better visible in a larger u- window in the inset. 
Such negative curvature is at odds with the w-dependence 
of AE'^si coming from the contribution alone of as 
predicted in Eq. ([?7)) . This suggests that the negative 
contribution from 71 is here of the same order of the term 



A 



On the other hand, the ai term alone would 



give an extrapolation of the ratio /S.E/ui at — > of 
the order Yuriu^^o /S.E / u^^ « —0.73 eV/A^ much bigger 
than what observed As a matter of fact, we can nicely 
reproduce the DFT data by taking = 0.30 eV/A. The 
comparison between DFT calculations and the TB model 
with this value of 04 reasonably good, as shown in Fig. 
[71 We have to stress however that, unlikely the other 
parameters ai that were obtained in a direct unbiased 
way by a high-precision extrapolation for u 0, since 
a4 was deducted in an indirect way from the knowledge 
of ai, and given the numerically scattered DFT data in 
Fig. [71 this value = 0.30 eV/A must be considered 
just as an indicative electron-phonon coupling for this 
hopping parameter. 



In this paper we have employed a combined TB and 
DFT approach to evaluate the deformation potential 
in single-layer and bilayer graphene associated with the 
modulation of the different hopping parameters. In or- 
der to avoid any fitting procedure, we have focused on the 
low-energy levels ty at the high-symmetry point K and 
we have characterized the electron-phonon coupling ai 
for each hopping term by a careful analysis of the frozen- 
phonon dependence of upon the lattice displacement 
for different lattice modes. In this way we were able to 
determine within a unique framework all the deforma- 
tion potentials ai for both the intralayer (i = 0) and 
interlayer hoppings (i = 1, 3, 4) as well as the TB param- 
eters 71, S. We summarize in Tabl^Jour results for a^. 
We can also compare these values with the estimates of 
the absolute magnitude of the corresponding hopping pa- 
rameters, as reported in the right column in Table [H The 
correlation between these two quantities is also shown in 
Fig. [21 which reveals an almost perfect linear scaling of 
ai with 7i. A mean-square fitting procedure gives 



(28) 



where A = 0.141 eV/A and B = 1.365 A^i. We would 
like to stress the importance of such robust underlying 
correlation between the magnitude of the hopping term 
and the corresponding electron-phonon interaction inde- 
pendently on the precise value of 7^. It is indeed well 
known that the estimates of the hopping parameters 7^ 
can significantly depend on the fitting procedure as well 
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i 


m (eV/A) 


M (eV) 


(IL) 


4.41 


3.12* 


(2L) 


4.40 


3.12* 


1 (2L) 


0.61 


0.351^ 


3 (2L) 


0.54 


0.29* 


4(2L) 


0.30 


0.12* 



* From Ref.-*^ 
^ present work 



TABLE I: Electron-phonon coupling at associated with each 
hopping parameter 7^ in single layer (IL) and bilayer (2L) 
graphene. We also show, in the right column, representative 
values of the TB parameters 71. We provide an estimate of 
71, while ji for i = 0, 3, 4 are taken from Ref. [4^ . 
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tally the size of the electron-phonon coupling once the 
band parameters 7^ are extracted experimentally, for 
instance by means angle- resolved photoemission spec- 
troscopy (ARPES). In particular, the evolution of the 
electron-phonon coupling can be followed as a function 
of doping, applied electric-field, strain, etc... This can be 
can in a quite easy and safe way for 70, by looking at the 
linear conical dispersion at the K point, and for 71, by 
looking at the upper and lower band energy splitting at 
the same K point in bilayer graphene and graphite. Ex- 
perimental determinations of 73 and 74 have been also 
provided in literature. 

Our analysis provides thus a crucial, and previously 
missing, information to include quantitatively the role of 
the lattice deformations on the electronic, transport and 
optical properties of multilayered graphene. The effects 
of the lattice deformations on the electronic structure 
can be included in TB models involving the deformation 
potential associated with higher hopping terms than the 
nearest-neighbor ones. 
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Appendix A: Dirac cone reconstruction upon E2g, 
lattice distortion 



FIG. 8: Plot of the ai vs. 7^ parameters obtained from dif- 
ferent approaches. Empty circles are data obtained by the 
present work where Oi was estimated by the frozen-phonon 
technique and 7i, when not available, were taken from Ref. 14^ . 
Filled squares are data collected by Ref. [3^ using a wide vari- 
ety of techniques, including Hartree-Fock, the hybrid B2LYP 
functional, LDA, GGA and GW. Inset; same data on a larger 
scale. 



as on the inclusion of many-body effects in first-principles 
band structure for example, within the GW scheme. A 
detailed study of this issue, including also Hartree-Fock 
(HF) calculations, is provided in Ref. I38ll45l . where they 
also estimate within the same level of approximation the 
overall electronic 7r-bandwidth, related to 70, and the 
electron-phonon coupling in single-layer graphene and 
graphite. Their results are also plotted in Fig. [51 where 
we have translated the high-energy 7r-band splitting Aeu 
at the M point in the hopping parameter through the phe- 
nomenological relation Acm = I-2l7o. Also in this case, 
considering the widest variety of approaches (HF, LDA, 
GGA, hybrid B3LYP functional and GW), the trend is 
almost perfectly linear. 

Apart the fundamental implications of this result, it 
suggests a well, defined way to estimate experimen- 



In this Appendix we discuss in more details the origin 
and the phenomenology of the reconstruction of the Dirac 
point at the K edge for a critical value of the i?2g2 lattice 
distortion, as pointed out by DFT calculations in Fig. [S] 
and confirmed by the TB model. 

As a starting point we remind that in realistic undis- 
torted bilayer graphenes, electronic processes like the 
"skew" hopping 73 split the the parabolic Dirac cone 
in four linear Dirac points4^ In the simplest TB model 
with only 70-71-73 hoppings, the four Dirac points are lo- 
cated respectively at k = (0,0), (^3,0), (— A:3/2, •\/3A;3/2), 
(-fc3/2,-\/3fc3/2), where fcs = 71 73/70 ^iuf-— 

In order to investigate the role of the E2g2 lattice dis- 
tortion, we expand the Hamiltonian (jlip for small but 
finite k = (k^jky). Neglecting here for simplicity the 
terms 74, 5 that break the particle-hole symmetry, we 
can thus write: 











^3^lu 







71 








71 





TTO.u 







^0,u 






(M) 



where 7ro,„ = kx + iky + OQUx, tts^u = kx+iky-asu^, and 
where aq = Sao/hvp, 03 = Sa^cosO/hvp, 71 = ji/hvp- 
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The evolution of the Dirac points, corresponding to 
the low-energy states of (jA2| , as a function of Ux , in the 
relevant region < 0, is shown in Fig. [9l The inset 
shows also a zoom close to the K point. In similar way 
as it has been reported for uniaxial strain,^ also upon 
the optical i?2g2 lattice distortion the Dirac points shift 
away from their original location for u = 0. While such 
shift is monotonic for the three "leg parts", the shift of 
the central one is however non monotonic, with a initial 
departure from the K point, followed by a turn back along 
the opposite direction. Hence, at a critical value = 
—7103/00 the "central part" will eventually cross again 
the K point and then continue moving on the opposite 
side. 

We can quantify this evolution by focusing on the axis 
kx and tracing the evolution of the roots of Eq. (|A2[) for 
ky = 0. A straightforward analysis gives thus: 



FIG. 9: (color online) Evolution of the low-energy Dirac-like 
states (e « 0) as a function of the -E292 lattice displacements 
Ux- Inset: a zoom in the region \kJ, \ky \ < 0.001 A-\ Colors 
refer here to the energy distance from the Dirac points at 
e = 0, The color scale has been adapted in each panel to 
make more visible the low energy states, with e = being the 
darked regions. 



In the absence of particle-hole asymmetry, the four 
Dirac cones lie at the same energy e = also in the 
presence of lattice distortion. We can thus trace their 
evolution as a function of by analyzing the solution 



det 



(Ux) 



0. 



(A2) 



kxA 



71 W3 - 2aQUx 



±^\J^lvl - 47iW3aoUx - 47ia3M:r, (A3) 



where kx^- is the non monotonic solution for Uj; < 
starting from k — (0, 0) at = and kx^+ is the sec- 
ond shifting away solution starting from k = (^3,0). 
From Eq. (|A3|) we thus get a critical value Ux = 
—7103/00 — — a3cos6'7i/3ao. Similar calculations can 
be generalized including the crystal field S which breaks 
the particle-hole symmetry. We get in this case the result 
Ux = —as cos6'(7i ± S)/3aQ, as reported in Sec. Ill B 21 
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